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ABSTRACT 

Multidimensional reactive flow models of accreted hydrogen rich envelopes 
on top of degenerate cold white dwarfs are very effective tools for the study of 
critical, non spherically symmetric, behaviors during the early stages of nova 
outbursts. Such models can shed light both on the mechanism responsible for 
the heavy element enrichment observed to characterize nova envelope matter and 
on the role of perturbations during the early stages of ignition of the runaway. 
The complexity of convective reactive flow in multi-dimensions makes the 
computational model itself complex and sensitive to the details of the numerics. 
In this study, we demonstrate that the imposed outer boundary condition can 
have a dramatic effect on the solution. Several commonly used choices for the 
outer boundary conditions are examined. It is shown that the solutions obtained 
from Lagrangian simulations, where the envelope is allowed to expand and 
mass is being conserved, are consistent with spherically symmetric solutions. In 
Eulerian schemes which utilize an outer boundary condition of free outflow, the 
outburst can be artificially quenched. 



Subject headings: methods: hydrodynamics , numerical , (stars:) novae 



1. Introduction 

The standard model for nova outbursts is a thermonuclear runaway in an accreted 
hydrogen envelope on the surface of a white dwarf. On time scales of days prior to the 
runaway, the flow in the partially degenerate hydrogen rich envelope is convectively unstable. 
The convective instability is driven by nuclear burning that is highly sensitive both to the 
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temperature and to the chemical abundances. The details of the abundance mixing within 
and at the edge of the convective zone and the fate of temperature fluctuations are therefore 
crucial elements of the runaway mechanism. In order to capture these effects, several 
multidimensional simulations of nova outbursts have been performed (Shankar, Arnett & 
Fryxell 1992, Shankar & Arnett 1994, Glasner & Livne 1995, Glasner, Livne & Truran 1997, 
Kercek, Hillebrandt & Truran 1998, Kercek et al. 1999, Dursi et. al. 2002,Alexakis et. al. 
2004). 

In this study we attempt to validate multi dimensional schemes by discriminating 
numerical effects from physical ones in the explosion phase. One problem is to identify 
proper initial conditions for the simulations. Since the pre-runaway quasi-hydrostatic 
accretion time is orders of magnitude longer than the runaway time, all of the published 
multi dimensional nova calculations start from a spherically symmetric ID model, mapped 
to 2D (or 3D). The computed measurable variables are sensitive to the adaptation of the 
ID initial model to the multidimensional model, a problem that was discussed by (Kercek, 
Hillebrandt & Truran 1998). 

Another problem is the implementation of a correct outer boundary condition. Nariai, 
Nomoto & Sugimoto (1980) discuss the sensitivity of ID Lagrangian models to the outer 
boundary conditions during the expansion phase. They claim that fine zones are needed 
in order to resolve the details of the expansion of the nova atmosphere from white dwarf 
dimensions to red giant dimensions. 

Two independent studies - those of (Glasner & Livne 1995, Glasner, Livne & Truran 
1997) and (Kercek, Hillebrandt & Truran 1998; Kercek et al. 1999) - used the same ID 
initial models but were led to different conclusions about the strength of the runaway and 
its ability to reproduce a fast nova. In order to identify the reasons for these differences, we 
carefully analyzed the sensitivities mentioned above. We argue here that the early stages of 
the flash phase, for which the evolution is still almost quasi-static, are highly sensitive to 
the outer boundary conditions. The sensitivity addressed in this study takes place during 
the explosion (flash) phase and the main issue is the difference between multi-dimensional 
numerical schemes (Lagrangian vs. Eulerian), which has no ID counterpart. The sensitivity 
is certainly not a problem of resolution (see next section). In the explosion (flash), phase 
the issue of resolution is important mainly in the mixing layers, but not outside at the 
atmosphere. In this phase, as the flash develops, the radial dimension of the envelope can 
increase by a factor of three (Fig. 3), affecting the pressure at the bottom of the accreted 
envelope where burning takes place. 

The question of proper boundary conditions for multi dimensional simulations during 
this phase is therefore, the main topic addressed in this work. We show that different, 
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commonly used, mass outflow conditions at the outer boundary can make the difference 
between a successful runaway model and a model for which the runaway is artificially 
quenched. We analyze the cause for this failure and conclude with a recommendation for 
boundary conditions that are consistent with the physics of the problem. 

2. Problematic outer boundary conditions - which and why 

Our initial ID model is a 1.14 M CO WD that has accreted a 2.6 x 1(T 5 M envelope 
(solar abundance) up to the stage at which a TNR develops. The model is composed of 130 
CO core zones and 50 accreted hydrogen zones. In ID, the main features of the runaway, 
such as the accreted mass up to the runaway, the peak temperature and the peak energy 
generation rate, are converged at that spatial resolution. Increasing the number of zones in 
the accreted hydrogen layer to 250 change those values by a fraction of a percent. 

When the temperature at the base of the envelope has reached 1 x 10 8 K (just before 
the ID runaway starts), the outer 6 x 10 -5 M of the star is mapped, without any change 
in the zoning, as a multidimensional model (2D or 3D) onto the multidimensional hydro 
solver. The evolution in 2D is followed without any initial perturbation. The ID code and 
the 2D code use the same difference scheme for the momentum equation, the same radial 
zoning and the same EOS. Therefore, after the initial mapping the envelope is almost 
everywhere in full hydrostatic equilibrium. Deviations from hydrostatic equilibrium at local 
points are less then one part in ten thousand. For models at earlier stages in the runaway 
where the envelope is still stable against convection, the 2D mapped model stays in perfect 
hydrostatic equilibrium for a hundred seconds (the typical time step, dictated by the sound 
crossing time is about 10~ 3 seconds). In this way we avoid the problems of mapping the ID 
models to multi dimensional models, a problem that is thoroughly discussed by (Zingale et. 
al. 2002). 

The numerical scheme of the ID model is a Lagrangian scheme for which mass at the 
outer boundary is conserved and there is no overshoot mixing between the core matter and 
the accreted envelope. The multidimensional models try to analyze closely the details of 
the convective flow and are therefore basically Eulerian. Because of the extreme effect that 
overshoot mixing at the inner interface has on the nuclear burning in 2D it is impossible to 
judge the accuracy of the 2D model by the ID results. 

The models presented here are all 2D models. The 2D hydro VULCAN solver (Livne 
1993 consists of two computational stages within one time step. As a first stage a purely 
Lagrangian hydro time step is performed (explicit in our case). In the second stage , a 
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high-order explicit remapping is performed. In this remapping stage one can remap the new 
flow field onto any desired grid; the original Eulerian grid, or any other adaptive grid. The 
scheme as a whole is therefore called hereafter Arbitrary Lagrangian Eulerian (ALE). The 
lateral boundary conditions at the sides of the multidimensional slice can be either fixed or 
periodic with only a slight effect on the solution. 

The computational slice includes a significant part of the static, inert, white dwarf 
core. The deep cold part of the core justifies a fixed inner boundary condition. We assumed 
initially that, since the computational slice covers a few pressure scale heights, there will 
not be any sensitivity to the exact outer boundary condition. Subsequently, we determined 
that this assumption is quite unjustified. It is well known that the runaway is a threshold 
phenomena, the pressure at the base of the envelope must be in excess of a certain value of 
a few time 10 19 erg/cm 3 (Fujimoto 1982, MacDonald 1983). Spurious mass loss from the 
computational grid at the outer boundary slightly decreases the pressure at the bottom of 
the burning zone. In the worst case, this mass loss totally quenches the runaway while in 
other cases it leads to a very fast non-physical turnoff of the runaway. Since the physical 
mechanism for the turnoff is the decrease in pressure due to expansion, it is clear that 
a correct boundary condition is essential. In order to demonstrate the problem and its 
solution, we present here the results we have obtained for four different schemes for the 
outer boundary condition that were tested: 

• Eulerian scheme with free outflow at the outer boundary (Eulerian - open in Fig. 1): 
For this scheme, matter with positive radial velocity in the outer zone is lost from the 
computational grid at each step. 

• Eulerian scheme with modified mass flux outer boundary condition (Eulerian 
inflow+outflow in Fig. 1): This scheme derives from the observation that during the 
development of the runaway the outer boundary is convective. It follows that at any 
time some of the outer grid points have velocities directed outward and some have 
incoming velocities. Therefore, as long as the evolution is quasi-static, the amount 
of mass that leaves the last original zone should flow back onto the grid once the 
convective turnover of this convection cell is completed. Based on this assumption, 
we added a modified mass flux outer boundary condition. For this scheme, whenever 
the velocities in the outer zones point outward, the mass flux is calculated and mass 
is lost. When the velocities point inward, however, we impose an incoming mass flux 
with the assumption that the outer zone is also the donor. This method is in a way 
equivalent to models working with "ghost cells" for extrapolation of the mass fluxes 
in the case of "open boundaries" . 



- 5 - 



• Eulerian scheme with no mass flux (fixed) outer boundary condition (Eulerian fix 
- closed in Fig. 1): For this option, matter is forced to stay within the grid (no 
mass flux is allowed to pass into or out of the outer edge of the grid). Mass that is 
artificially accumulated near the outer boundary produces extra outer pressure. This 
extra pressure is only an artifact of the numerical model. 

• An Arbitrary Lagrangian Eulerian (ALE) scheme for which the radius of the outer 
zone is defined in a way that conserves mass (Lagrangian in Fig. 1). This scheme 
enables us to overcome three major obstacles that the nova models face: a) The 
Lagrangian hydro step takes care of the dimensional expansion of the computed 
regime, b) The remap stage enables us to keep a radial grid all the time. The radial 
grid has K_max zones. The number of radial zones is fixed and they are numbered 
by the index k, k = 1, Kjmax . The transversal zones are numbered by the index 
j, j — 1, Jjmax. All the transversal zones (j = 1, Jjnax) for each k are defined to 
be the k-th shell, c) The arbitrariness of the Eulerian grid enables us to model the 
burning zones at the bottom of the hydrogen rich envelope with very delicate zones in 
spite of the dimensional expansion due to the runaway. 

In order to realize these goals at each time step, a "radial" mass conserving grid is 
defined ( after the Lagrangian part of the step) in the following way. For each k shell 
of grid cells that started the step with the same radius a new radius is defined so 
that the mass of that shell is conserved. In the new grid, defined for the next step, 
the lateral boundaries of each grid cell are kept fixed (Jjnax lateral equal zones). 
The radial coordinate of the whole k-th shell of grid cells is defined by interpolation 
between the initial Eulerian radial grid and the new mass conserving grid. The only 
restriction is that, at the outer parts of the grid, the K_max shell will follow the new 
mass conserving grid so that the total mass is conserved. Since we usually want to 
keep the fine zoning in the inner parts of the grid where intense burning and mixing 
takes place ( item (c) above), the new grid smoothly interpolates between the initial 
Eulerian grid in the inner zones and the "Lagrangian" grid at the outer zones. In this 
way the outer boundary is expanding with time. The interpolation we use is of the 
following form: 

R new (k) = f x R Lag (k) + (1 - /) X R (k); f = (k/K m ax) n (1) 

where: 

— A; is the index of the radial k-th shell. 

— RLag(k) is the radius of the k-th shell in the mass conserving Lagrangian grid. 
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- Ro(k) is the original radius of the k-th shell in the initial Eulerian grid 

— n is a power low. The default value we use is n = 3. 

In order to avoid small differences that might appear during the initial stages, (just 
after the mapping of the ID model into 2D), all four models started from a single 2D profile 
of the ALE scheme 20 seconds after it was mapped into 2D. The convective cells are already 
fully developed at this stage. 

With the present resolution, it is difficult to obtain any significant information on the 
light curve. The results concerning the strength of the outburst are therefore judged only 
by comparison with ID spherically symmetric models (Nariai, Nomoto & Sugimoto (1980) 
; Prialnik & Kovetz 1984 ; Starrfield, Sparks & Truran 1985 ; Starrfield, Sparks & Truran 
1986 ; Prialnik 1986). In this respect, we define a successful nova to be a TNR for which 
the amount of thermonuclear energy pumped into the envelope on a dynamic time scale is 
comparable to the binding energy of the envelope. For the case we have considered of a 
1.14 M WD with an accreted envelope mass of 2.6 x 10~ 5 M© , the binding energy of the 
envelope is approximately 10 46 erg. The criterion we therefore adopt for a successful nova 
is an overall reaction rate greater then 1.0 x 10 44 erg/sec for more then 50 seconds (typical 
dynamical time of the accreted envelope). We examine the sensitivity of the 2D evolution 
of the models by comparing the following measurable quantities; 

• The time history of the overall energy production rate (Fig. 1). 

• The time history of the pressure at the base of the computational slice (Fig. 2). 

The ALE scheme is expected to be the most accurate. Indeed, for this scheme 
the energy production rate increases continuously and a runaway is taking place. The 
characteristics of this model throughout its evolution is compatible with the ID evolution of 
CO enriched models. We have therefore chosen to compare all other models to this model. 

The free Eulerian-open scheme gives a significantly different result: essentially no 
runaway. Although the initial pressure at the outer zone is two orders of magnitude less 
than the pressure at the bottom zone, the small decrease in pressure on a few dynamic times 
(the flow is extremely sub-sonic) quenches the runaway. The differences in the pressure 
profiles are evident in Fig. 2. In Fig. 3 we demonstrate the major difference between the 
two schemes. The figure shows a color map of the logarithm of the energy production rate 
within the envelope close to the maximum of energy production (time=100 sec, see Fig. 1) 
in both models. For the ALE Lagrangian mass conserving scheme we see the high energy 
production rates at the bottom, reaching values of 1 — 2 x 10 17 erg/gr/sec. For the pure 
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Eulerian scheme the rates are below 1 x 10 15 erg/gr/sec. The differences in the geometric 
dimensions of the envelope in both models is clear, once observed it confirms the necessity 
of the Lagrangian approach. 

The modified Eulerian scheme conserves mass for a while, as long as the matter is 
convecting in the absence of overall expansion. Once the runaway starts to take place (at 
about t = 80 seconds), however, the net expansion velocities in the outer regions are higher 
than any component of the convective flow at this zone. Subsequently, mass loss quenches 
the runaway (Fig. 1). 

The Eulerian-closed scheme agrees quite well with the ALE scheme until the point at 
which pressure starts to drop in the Lagrangian scheme (t = 180 seconds in Fig. 2). The 
agreement is both from the point of view of energy production rates (Fig. 1) and of the 
topology of the convective cells. Once we prevent the expansion of matter at the outer 
boundary, however, it flows back and artificially enhances the convective velocities and the 
mixing (Fig. 2). This scheme thus extends the runaway artificially so that in the absence 
of expansion there will be no turnoff of the TNR. 

3. Discussion and Conclusion 

The onset of the thermonuclear runaway in nova outbursts occurs once the accreted 
envelope reaches a critical pressure. A detailed knowledge of the evolution of the pressure 
at the base of the envelope is therefore essential to our understanding of the development of 
the runaway. The numerical simulations presented in this paper show the major role which 
the outer boundary conditions play in the evolution of the outburst. Of the four schemes 
we examine in this paper, only the ALE scheme provides a correct physical solution to the 
problem. The free outflow (open) boundary condition used in Eulerian schemes artificially 
quenches the runaway. Use of such a boundary condition in Eulerian simulations can explain 
the main differences between the models presented by Glasner, Livne & Truran 1997 and 
those of Kercek, Hillebrandt & Truran 1998. The fact that in the latter simulation virtually 
all of the hydrogen disappears from the grid at late times (Fig. 8 in Kercek, Hillebrandt & 
Truran 1998) supports this assumption (Fig.3). 

For the same reason, the modified boundary condition we have described fails to 
reproduce the physical features of nova outbursts during the expansion phase. It can serve 
properly to represent the boundary conditions only in the case of a quasi static convective 
envelope. 

A possible remedy to this problem could be an addition of empty grid cells at the outer 
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part of the computed region ("ghost cells"). The amount of expansion demonstrated in 
Fig. 3 shows that this solution is not appropriate for the problem presented here. 

To conclude, we state that proper multi dimensional simulations of Novae-like explosions 
should maintain the Lagrangian nature of the expanding envelope. The sensitivity we 
encountered is expected whenever the onset of a runaway is critically dependent upon the 
evolution of the pressure in a subsonic regime. 



Acknowledgments 

This work is supported in part at the University of Chicago (JWT) by the U.S. 
Department of Energy, under Grant B523820 to the ASCI Alliances Center for Astrophysical 
Flashes and Grant DE-FG02-91ER40606 in Nuclear Physics and Astrophysics and by 
the NSF under grant PHY 02-16783 for the Physics Frontier Center "Joint Institute for 
Nuclear Astrophysics." Ami Glasner, wants to thank the Department of Astronomy and 
Astrophysics, University of Chicago and the members of the ASCI Center for Astrophysical 
Thermonuclear Flashes for their hospitality during his visits to Chicago, where part of this 
work was done. 



- 9- 



REFERENCES 

Alexakis,A. , Calder,A.C. , Heger,A. , Dursi,L.J. , Truran,J.W. , Rosner,R. , Lamb,D.Q. , 
Timmes,F.X. , Fryxell,B. , Zingale,M. , Ricker,P. , and 01son,K. , ApJ,602,931. 

Dursi,L.J., Calder,A.C. , Alexakis,A. , Truran,J.W. , Zingale,M. , Fryxell,B. , Ricker,P. 
, Timmes,F.X. and 01so,K. 2002 in International Conference on Classical Nova 
Explosions. Sitges, Spain 2002, eds. Hernanz,M. & Jose, J. (AIP), p. 139. 

Fujimoto,M.Y.,1982,ApJ, 257, 752. 

Glasner,S.A. & Livne,E. 1995, ApJ, 445,L149 

Glasner,S.A., Livne,E., & Truran,J.W., 1997, ApJ,475,754. 

Kercek,A., Hillebrandt,W. & Truran,J.W., 1998, A&A, 337, 379 

Kercek Hillebrandt & Truran. 1999, A&A, 345, 831 

Livne,E., 1993, ApJ, 412, 634 

MacDonald,J.,1983,ApJ, 267, 732. 

Nariai,K., Nomoto,K., & Sugimoto,D. , 1980, PASJ, 92,473 
Prialnik, D. & Kovetz,A., 1984 ApJ, 281,367 
Prialnik, D., 1986 ApJ, 310,222 

Shankar,A., Arnett,D., & Fryxell,B.A., 1992, ApJ, 394,L13 
Shankar,A. & Arnett,D., 1994, ApJ, 433,216 
Starrfield,S. Sparks,W.M. & Truran,J.W., 1985, ApJ, 291,136 
Starrfield,S. Sparks,W.M. & Truran,J.W., 1986, ApJ, 303,L5 

Zingale,M. , Dursi,L.J. , ZuHone,J. , Calder,A.C. , Fryxell,B. , Plewa,T. , Truran,J.W. , 
Caceres,A. , Olson, K. , Ricker,P.M. , RileyK. , Rosner,R. , Siegel,A. , Timmes,F.X. 
, Valdimirova,N. , ApJS, 143,539. 



This preprint was prepared with the AAS IATeX macros v4.0. 



-10- 



Fig. 1. — Total burning rate as a function of time for all models (see text) 

Fig. 2. — Pressure at the base of the envelope as a function of time for all models (see text) 

Fig. 3. — Color map of the thermonuclear energy production rate at time=100 sec for the 
pure Eulerian case-right and the ALE Lagrangian scheme- left. The spatial coordinate is in 
units of 100 Km. The energy production rate is in erg/gr/sec. The rate scale is different in 
the two cases (see scale to the right of each model) (see text) 
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